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We present a computational method to accurately calculate Raman spectra from first principles 
with an at least one order of magnitude higher efficiency. This scheme thus allows to routinely 
calculate finite-temperature Raman spectra “on-the-fly” by means of ab-initio molecular dynamics 
simulations. To demonstrate the predictive power of this approach we investigate the effect of 
hydrophobic and hydrophilic solutes in water solution on the infrared and Raman spectra. 


I. INTRODUCTION 

Ah initio simulations of vibrational spectra, where the 
electronic degrees of freedom are explicitly taken into ac¬ 
count, often provides important insights into the struc¬ 
ture and dynamics of complex systems. Therefore, com¬ 
puter simulations nowadays represent an invaluable tool 
to rationalize and complement experimental measure¬ 
ments. The key quantity to compute the vibrational 
spectrum of a system is its dipole moment, which is the 
tendency of charge distribution towards inhomogeneity. 
From this point of view. Maximally Localized Wannier 
Functions (MLWFs) are particularly useful, since they 
allow to partition the total electronic density into indi¬ 
vidual fragment contributionsThese Wannier functions 
are defined as 
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where '0rnk(i*) are Bloch functions, R is a Bravais lat¬ 
tice vector, and V is the real-space primitive cell volume, 
while the integral is computed over the whole Brillouin 
zon^. The unitary J x J matrix Umn is periodic with 
respect to the wave-vector k, while '0mk(i*) are the eigen¬ 
states of a system as obtained by an electronic structure 
method, such as density functional theory (DFT)I^. To 
compute MLWFs, the total spread functional 

S = ^ ^ Sn —'^n | '^n) ^ : (2) 
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is minimized by appropriately chosen unitary rotations 
^4^2Thereof, one can use the expectation value of the 
periodic position operator r in the Wannier representa¬ 
tion, in order to find the centers of the localized func¬ 
tions for arbitrary symmetriesAs such, the polariza¬ 
tion of the electronic charge in a crystal, which in general 
has a periodic continuous distribution, can be unambigu¬ 
ously partitioned into localized contributionsHaving 
the centers and the associated spreads of the MLWFs and 


considering each MLWF as a charge distribution in space, 
it is not only possible to calculate the molecular and/or 
total dipole moments of a system, but also the corre¬ 
sponding time-correlation functions by means of ah-initio 
molecular dynamics (AIMD) simulations. The temporal 
Fourier transform of the autocorrelation between the to¬ 
tal dipole moments M(t) at different times, t, is propor¬ 
tional to the infrared (IR) absorptivity, which after 

e mploy ing the harmonic approximation, can be expressed 
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where p and n{v) are the frequency and the index of re¬ 
fraction, respectively, while denotes the ensemble- 

average in classical statistical mechanics. By applying a 
periodic electr ic field,^!^®^ usually using the Berry phase 
approachit is possible to obtain the polarizability 
tensor A via 


Aij — 


dM,(E) 

dEj 
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where Ej denotes the Cartesian component j of the ap¬ 
plied electric field, while Mi is the Rh component of the 
total dipole moment. However, calculating the deriva¬ 
tives numerically based on density functional perturba¬ 
tion theory (DFPT)[43tini or using higher-order finite dif¬ 
ference (FD) method^^is computationally rather expen¬ 
sive. The mean polarizability A = 1/3 Tr[A] is the quan¬ 
tity that is usually measured in experiment. With A 
known, one can obtain the isotropic Raman spectrum 
through the calcu lation of the autocorrelation between 

the polarizabilitie^^illiJ 
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In this paper we present a novel computational technique 
to efficiently calculate Raman spectra, where the polariz¬ 
ability of a system is represented as a sum over Wannier 
function polarizabilities denoted as a function of Wannier 
function volumes, and consequently by their spread. 
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II. WANNIER POLARIZABILITY METHOD 

Specifically, it has been shown that the molecular po¬ 
larizability change linearly with the volume of the elec¬ 
tronic cloud around a moleculeP^ Due to the fact that 
the electronic properties play the main role in quantifying 
the polarizability, it is reasonable to assume that the total 
isotropic polarizability of the system can be obtained as 
a sum over the Wannier function polarizabilities. The es¬ 
sentially same idea has been recently used to combine the 
quantum harmonic oscillator model with Wannier func¬ 
tions to calculate the non-local dynamic electron correla¬ 
tion due to oscillating charges assigned to MLWFs.l^ As 
such, the polarizability assigned to ith MLWF is given 
by 

= /3Sl ( 6 ) 

where Si is the spread of the ith MLWF and /3 is a pro¬ 
portionality constant. Therefore, the mean polarizability 
of the system can be written as 
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At variance to the FD technique, the numerical effort re¬ 
duces from six to just one single-point calculation plus an 
additional MLWF computation. This is to say that our 
novel MLWF-based method, which hereafter we will refer 
to as Wannier polarizability (WP) method, is at least five 
times more efficient than the conventional FD app roach. 
The speed-up with respect to DFPT is similar.^ 

The parameter p is determined by minimizing the 
Mean Absolute Relative Error (MARE) of the mean po¬ 
larizabilities with respect to reference calculations us¬ 
ing the ED approach. All of our DET calculations 
were conducted using the CP2K/Quickstep cod^^ in 
conjunction with a very accurate TZV2PX Gaussian 
basis sefP^ and the Perdew-Burke-Ernzerhof exchange- 
correlation functionaP^plus a damped interaction poten¬ 
tial to approximately account for long-range dispersion 
interactions.!^ The mean polarizabilities per molecule 
for various water clusters containing 2 to 13 molecules 
as obtained using the ED scheme are shown in Eig. 
The structures of the water clusters at Hartree-Eock 
level of theory were taken from the Cambridge Cluster 
Database.!^ The obtained results are in very good agree¬ 
ment with previously calculated pol arizab ilities for the 
same systems at DET level of theory.!^^^^ 

The optimized value of (3 for our WP method at min¬ 
imum MARE is p = 0.90. The corresponding isotropic 
polarizabilities of the water clusters are also displayed in 
Fig-B Even though the absolute polarizabilities of the 
ED and WP methods slightly differ for the smallest water 
clusters, the qualitative behavior with respect to cluster 
size is similar, which immediately suggest that the de¬ 
duced Raman spectra most likely differ only in their ab¬ 
solute intensities. To assess the latter, we have compared 



FIG. 1. Mean polarizabilities per water molecule for various 
water clusters containing 2 to 13 molecules as obtained using 
the FD (black circles) and WP (red squares) methods. The 
error bars for the red squares are the absolute relative errors 
at each point. 




FIG. 2. Atomic configuration of cyclohexane (a) and cyclo- 
hexanedodecol (b). Cyan, red and white spheres denote car¬ 
bon, oxygen and hydrogen atoms, respectively. 


the Raman spectra of a single cyclohexane molecule in 
the gas phase, shown in Fig. [^a), using both FD and WP 
methods. The spectra were obtained “on-the-fly” from 
simulations using the second-generation Gar-Parrinello 
method, where both the density matrix and Umn{t) were 
propagated tog ether with the nuclei to further speed-up 
the calculationsHowever, for each AIMD step, the 
Wannier functions were re-localized using the scheme of 
Berghold et al. to obtain genuine MLWFs.!^ Nonetheless, 
together with the WF method, this results in a combined 
acceleration for the calculation of the Raman spectra of 
at least an order of magnitude. At first the system was 
equilibrated in the canonical ensemble at 300 K for 10 ps 
using a discretized time-step of 0.5 fs, before sampling the 
polarizabilities in the micro-canonical ensemble for addi¬ 
tional 10 ps. In the case of FD approach, at each time- 
step we applied an external electric field of 0.0001 a.u. 
intensity along the ±x, 3zy and iLz directions, to cal¬ 
culate the polarizabilities. The eventual Raman spectra 
are depicted in Fig. As can be seen, the agreement 
between the FD and WP methods is excellent, and no 
frequency shift can be observed. The peaks around 800 
cm“^ are typically attributed to G-G stretching and GH 2 
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FIG. 3. Isotropic Raman spectra of cyclohexane molecule in 
the gas phase, as obtained by the FD (red dashed line) and 
WP (black straight line) methods. 


rocking vibrations, while the peaks around 3000 cni“^ 
are usually assigned to symmetric and asymmetric CH 2 
stretches. In the inset of Fig. [^the Raman activity for 
frequencies less than 1500 cm“^ are shown, where the 
higher bands are due to other CH bending modes, while 
the lower ones originates from the torsion and deforma¬ 
tion of the carbon rin^^^MH^ Qur simulated spectra are in 
good agreement with other theoretical and experimental 
results. 


III. RESULTS AND DISCUSSION 


We demonstrate our novel WP method by investigat¬ 
ing the effect of hydrophobic and hydrophilic molecules 
on the IR and Raman spectra of liquid water. To that 
extend we consider two systems in water solution: a 
single cyclohexane molecule (Fig. [^a)) to represent a 
hydrophobic solute, and the relatively similar, but hy¬ 
drophilic cyclohexanedodecol molecule (Fig.|^b)), where 
all hydrogen atoms are replaced by OH groups. Here¬ 
after we refer to these two systems as CW and COHW, 
respectively. Both, CW and COHW contains 128 light 
water molecules per unit cell, where the water density 
is set to the experimental value at ambient conditions. 
Contrary to our previous calculations, a smaller DZVP 
basis set has been employed. Again, both structures 
were first equilibrated in the canonical ensemble at room 
temperature for ^15 ps, before the dipole moments and 
polarizabilities were sampled “on-the-fly” for 10 ps in 
the micro-canonical ensemble. The simulated IR and 
Raman spectra of CW and COHW systems are shown 
in Fig. The IR activities for the frequencies below 
650 cm”, at ^1600 cm“^, as well as the weak peak at 
^ 2200 cm“^ are due to libration and hydrogen bond 
bending and stretching, 0-H bending, and a combina¬ 
tion of 0-H bending and libration modes of bulk water, 



FIG. 4. IR and Raman spectra of GW (a) and GOHW (b) sys¬ 
tems. Black and red curves show infrared and Raman spectra, 
respectively. 


respectively.E^®^ The peak at ^800 cm“^ in the IR and 
Raman spectra of the CW and COHW systems origi¬ 
nate most likely from C-C stretching modes of the car¬ 
bon ring.l^^HUl Moreover, IR activities in the frequency 
range 950-1100 cm“^ can be assigned to C-0 stretching 
modes of cyclohexanedodecol, while the peak at ^1350 
cm“^ seems to be due to its C-O-H bending mode.S^By 
comparison of FigQa) with FigQb) we attribute the 
IR-active peak at r\j 3000 cm ^ of the CW system to C-H 
stretching modes of the cyclohexane molecule. The re¬ 
maining frequenc ies ab ove above 3000 cm“^ are the 0-H 
stretching modesHowever, at variance to bulk wa¬ 
ter, we observe a generally larger splitting between the 
symmetric and asymmetric stretching modes, which im¬ 
mediately suggests that the asymmetry of the hydrogen- 
bond network is more pronounced due to the presence of 
the solute.E^J^ In the case of the hydrophobic cyclohex¬ 
ane molecule the effect is even more pronounced. More¬ 
over, for the CW system the so-called dangling 0-H bond 
peak at ^3650 cm“^ is more distinct and gives rise to a 
shoulder as can be seen in Fig.[^a).I^^Jf^ We believe that 
the latter is a consequence of fleetingly broken hydrogen- 
bonds of the distorted hydrogen-bond network which is 
spanned around the hydrophobic solute. 


IV. CONCLUSIONS 


In summary, we have presented a novel method that 
allows to efficiently calculate IR and in particular Ra¬ 
man spectra “on-the-fly” within AIMD simulations. To 
that extend we exploit the fact MLWFs, which are at 
the core of this new approach, can be utilized to parti¬ 
tion the charge distribution of the system into localized 
fragments. Therefore, the total isotropic polarizability 
can be calculated as a sum over the Wannier polarizabil¬ 
ities, which are assumed to be proportional to its volume 
and determined by its spread. Together with an exten¬ 
sion of the second-generation Car-Parrinello method to 
propagate Umn{t) along with the nuclei, followed my a 
re-localization to obtain genuine MLWFs, a speed-up of 
one order of magnitude has been observed. Using this 
approach, we calculated IR and Raman spectra for cyclo¬ 
hexane and cyclohexanedodecol solutes in ambient bulk 
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water. We found that the former hydrophobic solute give 
rise to a shoulder at around ^3650 cm“^, which is due 
to momentarily dangling 0-H bonds. In any case, we 


conclude by noting that this development facilitates to 
routinely calculate finite temperature spectra with only 
minimal extra computational cost. 
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